##install.packages("lme4")
library("lme4")

##### full data set
setwd("/Users/nicola/Documents/PISA/Our replication/")

mydata <- read.csv("Pisa data.csv")

###The code below is for Seaton's basic BFLPE model and for the BFLPE moderated by individual ability. We ran these both with and without weights, although only the code for with weights is reported here. For without weights, simple delete weights = mydata$NewWeight in the lmer call. ZSCMAT is the dependent variable Maths self-concept. ZPV1MATH is the first plausible value for individual ability - there are 5 PVs, Quad_ZPV1 is simple square of ZPV1MATH.

###It is not possible to retrieve standard errors for random effects from lmer. This is an intentional omission on the part of the developer Douglas Bates. He offers an alternative by calculating confidence intervals using mcmcsamp, however this too failed because he has not written the appropriate code to deal with certain cases (including ours).

### It may be necessary to extend the number of iterations allowed. Note that the Basic Model with weights experienced troubles with convergence (note different outputs on RCE, non-RCE machines!) also a problem for HISEI --> to extend max iterations (default 300), include control = list(maxIter = 500) in lmer call

#####################################################################################
#####################################################################################
######### BFLPE basic model: code to run multi-level models for each standardised plausible value (ZPV)

five.model.summary <- matrix(data = NA,16,5)

for (i in 0:4) {

##note column 60 of mydata is ZPV1MATH (ie Linear Ability) and column 88 is Quad_ZPV1 (ie Quadratic Ability)

lme4.out <- lmer( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + (1 + mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility| COUNTRY) + (1 + mydata[,60+i] + mydata[,88+i] | newSchoolID), data = mydata,  weights = mydata$NewWeight)

##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates, then random effects estimates for schools (Intercept, ZPV, QZPV) then for country (Intercept, ZPV, QZPV, SAA) and then the final residual

five.model.summary[,i+1] <- c(summary(lme4.out)@coefs[,1],summary(lme4.out)@coefs[,2],as.numeric(summary(lme4.out)@REmat[,3]))

}
print(five.model.summary)

BFLPE <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability", "L2 school intercept", "L2 Linear Ability", "L2 Quad Ability","L3 country intercept", "L3 Linear Ability", "L3 Quad Ability","L3 School Average Ability", "L1 individualintercept"),c("Estimates",round(rowMeans(five.model.summary)[1:4], digits = 4), round(rowMeans(five.model.summary)[9:16], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[5:8], digits = 4),"","","","","","","","")))


#####################################################################################
#####################################################################################
######### BFLPE moderated by Ability: multi-level regression

five.model.summary <- matrix(data = NA, 20, 5)

for (i in 0:4) {
	
lme4.out <- lmer( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + mydata[,60+i]*SchoolAverageAbility + mydata[,88+i]*SchoolAverageAbility + (1 + mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility| COUNTRY) + (1 + mydata[,60+i] + mydata[,88+i] | newSchoolID), data = mydata,  weights = mydata$NewWeight)

five.model.summary[,i+1] <- c(summary(lme4.out)@coefs[,1],summary(lme4.out)@coefs[,2],as.numeric(summary(lme4.out)@REmat[,3]))
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA, SAA*LA, SAA*QA), then standard errors of fixed estimates, then random effects estimates for schools (Intercept, ZPV, QZPV) then for country (Intercept, ZPV, QZPV, SAA) and then the final residual

}
print(five.model.summary)

BFLPE_Ability <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability","SAA*LA","SAA*QA", "L2 school intercept", "L2 Linear Ability", "L2 Quad Ability","L3 country intercept", "L3 Linear Ability", "L3 Quad Ability","L3 School Average Ability", "L1 individualintercept"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4), round(rowMeans(five.model.summary)[13:20], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4),"","","","","","","","")))




####### Below is the code for running the models reported in Seaton's table's 2 and 3. We simply ran out of time to run these models - the first one ran overnight with no conclusion. However, the code is included for interest.

#####################################################################################
#####################################################################################
######### BFLPE moderated by Highest Parental Occupation, HISEI: multi-level regression

five.model.summary <- matrix(data = NA, 22, 5)

for (i in 0:4) {
	
lme4.out <- lmer( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZHISEI + SAA_ZHISEI + (1 + mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZHISEI| COUNTRY) + (1 + mydata[,60+i] + mydata[,88+i] + ZHISEI| newSchoolID), data = mydata,  weights = mydata$NewWeight)

five.model.summary[,i+1] <- c(summary(lme4.out)@coefs[,1],summary(lme4.out)@coefs[,2],as.numeric(summary(lme4.out)@REmat[,3]))
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates, then random effects estimates for schools (Intercept, ZPV, QZPV) then for country (Intercept, ZPV, QZPV, SAA) and then the final residual

}
print(five.model.summary)

BFLPE_HISEI <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability","HISEI","SAA*HISEI", "L2 school intercept", "L2 Linear Ability", "L2 Quad Ability","L2 HISEI","L3 country intercept", "L3 Linear Ability", "L3 Quad Ability","L3 School Average Ability", "L3 HISEI", "L1 individual intercept"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4), round(rowMeans(five.model.summary)[13:22], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4),"","","","","","","","","","")))


#####################################################################################
#####################################################################################
######### BFLPE moderated by Highest Parental Education, HISCED: multi-level regression

five.model.summary <- matrix(data = NA, 22, 5)

for (i in 0:4) {
	
lme4.out <- lmer( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZHISCED + SAA_ZHISCED + (1 + mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZHISCED| COUNTRY) + (1 + mydata[,60+i] + mydata[,88+i] + ZHISCED| newSchoolID), data = mydata,  weights = mydata$NewWeight)

five.model.summary[,i+1] <- c(summary(lme4.out)@coefs[,1],summary(lme4.out)@coefs[,2],as.numeric(summary(lme4.out)@REmat[,3]))
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates, then random effects estimates for schools (Intercept, ZPV, QZPV) then for country (Intercept, ZPV, QZPV, SAA) and then the final residual

}
print(five.model.summary)

BFLPE_HISCED <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability","HISCED","SAA*HISCED", "L2 school intercept", "L2 Linear Ability", "L2 Quad Ability","L2 HISCED","L3 country intercept", "L3 Linear Ability", "L3 Quad Ability","L3 School Average Ability", "L3 HISCED", "L1 individual intercept"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4), round(rowMeans(five.model.summary)[13:22], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4),"","","","","","","","","","")))


#####################################################################################
#####################################################################################
######### BFLPE moderated by Home Educational Resources, HEDRES: multi-level regression

five.model.summary <- matrix(data = NA, 22, 5)

for (i in 0:4) {
	
lme4.out <- lmer( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZHEDRES + SAA_ZHEDRES + (1 + mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZHEDRES| COUNTRY) + (1 + mydata[,60+i] + mydata[,88+i] + ZHEDRES| newSchoolID), data = mydata,  weights = mydata$NewWeight)

five.model.summary[,i+1] <- c(summary(lme4.out)@coefs[,1],summary(lme4.out)@coefs[,2],as.numeric(summary(lme4.out)@REmat[,3]))
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates, then random effects estimates for schools (Intercept, ZPV, QZPV) then for country (Intercept, ZPV, QZPV, SAA) and then the final residual

}
print(five.model.summary)

BFLPE_HEDRES <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability","HEDRES","SAA*HEDRES", "L2 school intercept", "L2 Linear Ability", "L2 Quad Ability","L2 HEDRES","L3 country intercept", "L3 Linear Ability", "L3 Quad Ability","L3 School Average Ability", "L3 HEDRES", "L1 individual intercept"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4), round(rowMeans(five.model.summary)[13:22], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4),"","","","","","","","","","")))


#####################################################################################
#####################################################################################
######### BFLPE moderated by Cultural Possessions, CULTPOSS: multi-level regression

five.model.summary <- matrix(data = NA, 22, 5)

for (i in 0:4) {
	
lme4.out <- lmer( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZCULTPOSS + SAA_ZCULTPOSS + (1 + mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZCULTPOSS| COUNTRY) + (1 + mydata[,60+i] + mydata[,88+i] + ZCULTPOSS| newSchoolID), data = mydata,  weights = mydata$NewWeight)

five.model.summary[,i+1] <- c(summary(lme4.out)@coefs[,1],summary(lme4.out)@coefs[,2],as.numeric(summary(lme4.out)@REmat[,3]))
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates, then random effects estimates for schools (Intercept, ZPV, QZPV) then for country (Intercept, ZPV, QZPV, SAA) and then the final residual

}
print(five.model.summary)

BFLPE_CULTPOSS <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability","CULTPOSS","SAA*CULTPOSS", "L2 school intercept", "L2 Linear Ability", "L2 Quad Ability","L2 CULTPOSS","L3 country intercept", "L3 Linear Ability", "L3 Quad Ability","L3 School Average Ability", "L3 CULTPOSS", "L1 individual intercept"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4), round(rowMeans(five.model.summary)[13:22], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4),"","","","","","","","","","")))

#####################################################################################
#####################################################################################
######### BFLPE moderated by Extrinsic motivation, INSTMOT: multi-level regression

five.model.summary <- matrix(data = NA, 22, 5)

for (i in 0:4) {
	
lme4.out <- lmer( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZINSTMOT + SAA_ZINSTMOT + (1 + mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZINSTMOT| COUNTRY) + (1 + mydata[,60+i] + mydata[,88+i] + ZINSTMOT| newSchoolID), data = mydata,  weights = mydata$NewWeight)

five.model.summary[,i+1] <- c(summary(lme4.out)@coefs[,1],summary(lme4.out)@coefs[,2],as.numeric(summary(lme4.out)@REmat[,3]))
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates, then random effects estimates for schools (Intercept, ZPV, QZPV) then for country (Intercept, ZPV, QZPV, SAA) and then the final residual

}
print(five.model.summary)

BFLPE_INSTMOT <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability","INSTMOT","SAA*INSTMOT", "L2 school intercept", "L2 Linear Ability", "L2 Quad Ability","L2 INSTMOT","L3 country intercept", "L3 Linear Ability", "L3 Quad Ability","L3 School Average Ability", "L3 INSTMOT", "L1 individual intercept"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4), round(rowMeans(five.model.summary)[13:22], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4),"","","","","","","","","","")))

#####################################################################################
#####################################################################################
######### BFLPE moderated by Intrinsic motivation, INTMAT: multi-level regression

five.model.summary <- matrix(data = NA, 22, 5)

for (i in 0:4) {
	
lme4.out <- lmer( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZINTMAT + SAA_ZINTMAT + (1 + mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZINTMAT| COUNTRY) + (1 + mydata[,60+i] + mydata[,88+i] + ZINTMAT| newSchoolID), data = mydata,  weights = mydata$NewWeight)

five.model.summary[,i+1] <- c(summary(lme4.out)@coefs[,1],summary(lme4.out)@coefs[,2],as.numeric(summary(lme4.out)@REmat[,3]))
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates, then random effects estimates for schools (Intercept, ZPV, QZPV) then for country (Intercept, ZPV, QZPV, SAA) and then the final residual

}
print(five.model.summary)

BFLPE_INTMAT <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability","INTMAT","SAA*INTMAT", "L2 school intercept", "L2 Linear Ability", "L2 Quad Ability","L2 INTMAT","L3 country intercept", "L3 Linear Ability", "L3 Quad Ability","L3 School Average Ability", "L3 INTMAT", "L1 individual intercept"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4), round(rowMeans(five.model.summary)[13:22], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4),"","","","","","","","","","")))

#####################################################################################
#####################################################################################
######### BFLPE moderated by Maths Self-Efficacy, MATHEFF: multi-level regression

five.model.summary <- matrix(data = NA, 22, 5)

for (i in 0:4) {
	
lme4.out <- lmer( formula = ZSCMAT ~ mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZMATHEFF + SAA_ZMATHEFF + (1 + mydata[,60+i] + mydata[,88+i] + SchoolAverageAbility + ZMATHEFF| COUNTRY) + (1 + mydata[,60+i] + mydata[,88+i] + ZMATHEFF| newSchoolID), data = mydata,  weights = mydata$NewWeight)

five.model.summary[,i+1] <- c(summary(lme4.out)@coefs[,1],summary(lme4.out)@coefs[,2],as.numeric(summary(lme4.out)@REmat[,3]))
##gives fixed effects estimates (Intercept, ZPV, QZPV, SAA), then standard errors of fixed estimates, then random effects estimates for schools (Intercept, ZPV, QZPV) then for country (Intercept, ZPV, QZPV, SAA) and then the final residual

}
print(five.model.summary)

BFLPE_MATHEFF <- data.frame(cbind(c("","Intercept", "Linear ability", "Quadratic Ability", "School Average Ability","MATHEFF","SAA*MATHEFF", "L2 school intercept", "L2 Linear Ability", "L2 Quad Ability","L2 MATHEFF","L3 country intercept", "L3 Linear Ability", "L3 Quad Ability","L3 School Average Ability", "L3 MATHEFF", "L1 individual intercept"),c("Estimates",round(rowMeans(five.model.summary)[1:6], digits = 4), round(rowMeans(five.model.summary)[13:22], digits = 4)),c("Standard Errors",round(rowMeans(five.model.summary)[7:12], digits = 4),"","","","","","","","","","")))


########################################
### Results tables

BFLPE
BFLPE_Ability
BFLPE_HISEI
BFLPE_HISCED
BFLPE_HEDRES
BFLPE_CULTPOSS 
BFLPE_INSTMOT 
BFLPE_INTMAT
BFLPE_MATHEFF
